06/02/2018

Example Data - Golfing

## # A tibble: 19 x 5
##     dist successes tries         p    error_sd
##    <int>     <dbl> <dbl>     <dbl>       <dbl>
##  1     2      1346  1443 0.9327789 0.006591874
##  2     3       577   694 0.8314121 0.014211556
##  3     4       337   455 0.7406593 0.020546548
##  4     5       208   353 0.5892351 0.026185062
##  5     6       149   272 0.5477941 0.030178131
##  6     7       136   256 0.5312500 0.031188905
##  7     8       111   240 0.4625000 0.032183960
##  8     9        69   217 0.3179724 0.031613007
##  9    10        67   200 0.3350000 0.033374766
## 10    11        75   237 0.3164557 0.030211036
## 11    12        52   202 0.2574257 0.030762402
## 12    13        46   192 0.2395833 0.030803744
## 13    14        54   174 0.3103448 0.035072250
## 14    15        28   167 0.1676647 0.028907578
## 15    16        27   201 0.1343284 0.024052622
## 16    17        31   195 0.1589744 0.026184896
## 17    18        33   191 0.1727749 0.027354921
## 18    19        20   147 0.1360544 0.028277490
## 19    20        24   152 0.1578947 0.029576394

Visualise the data

Logistic Regression?

logistic_binom <- glm(p ~ dist,
                      weights = tries,
                      data = golf,
                      family = binomial(link = "logit"))
logistic_binom
## 
## Call:  glm(formula = p ~ dist, family = binomial(link = "logit"), data = golf, 
##     weights = tries)
## 
## Coefficients:
## (Intercept)         dist  
##      2.2312      -0.2557  
## 
## Degrees of Freedom: 18 Total (i.e. Null);  17 Residual
## Null Deviance:       2411 
## Residual Deviance: 255.3     AIC: 365.9

… in Stan (the easy version)

library(rstanarm)
stan_logistic_binom <- stan_glm(p ~ dist,
                                weights = tries,
                                data = golf,
                                family = binomial(link = "logit"))
## 
## SAMPLING FOR MODEL 'binomial' NOW (CHAIN 1).
## 
## Gradient evaluation took 9.8e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.98 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.110557 seconds (Warm-up)
##                0.114992 seconds (Sampling)
##                0.225549 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'binomial' NOW (CHAIN 2).
## 
## Gradient evaluation took 2e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.2 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.112554 seconds (Warm-up)
##                0.105586 seconds (Sampling)
##                0.21814 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'binomial' NOW (CHAIN 3).
## 
## Gradient evaluation took 2.2e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.22 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.100435 seconds (Warm-up)
##                0.107553 seconds (Sampling)
##                0.207988 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'binomial' NOW (CHAIN 4).
## 
## Gradient evaluation took 2.2e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.22 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.113862 seconds (Warm-up)
##                0.090601 seconds (Sampling)
##                0.204463 seconds (Total)

… in Stan

coefficients(logistic_binom)
## (Intercept)        dist 
##   2.2312106  -0.2556919
coefficients(stan_logistic_binom)
## (Intercept)        dist 
##   2.2342119  -0.2560822

Priors

prior_summary(stan_logistic_binom)
## Priors for model 'stan_logistic_binom' 
## ------
## Intercept (after predictors centered)
##  ~ normal(location = 0, scale = 10)
## 
## Coefficients
##  ~ normal(location = 0, scale = 2.5)
##      **adjusted scale = 0.44
## ------
## See help('prior_summary.stanreg') for more details
stan_logistic_binom2 <- stan_glm(p ~ dist,
                                weights = tries,
                                data = golf,
                                family = binomial(link = "logit"),
                                prior = laplace(location = 0, scale = 10),
                                prior_intercept = cauchy(location = 0, scale = 2.5))

Predictions

predicted_proportions <- predict(logistic_binom, type='response')
str(predicted_proportions)
##  Named num [1:19] 0.848 0.812 0.77 0.722 0.668 ...
##  - attr(*, "names")= chr [1:19] "1" "2" "3" "4" ...
stan_linepred <- posterior_linpred(stan_logistic_binom, transform = T)
str(stan_linepred)
##  num [1:4000, 1:19] 0.849 0.846 0.839 0.848 0.842 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ iterations: NULL
##   ..$           : chr [1:19] "1" "2" "3" "4" ...

Posterior Predictive

logist_binom_predict <- posterior_predict(stan_logistic_binom, draws = 500)
str(logist_binom_predict)
##  'ppd' int [1:500, 1:19] 1231 1220 1196 1224 1245 1194 1261 1226 1224 1226 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:19] "1" "2" "3" "4" ...
# convert sucesses to proportions
logist_binom_predict <- t(t(logist_binom_predict) / golf$tries)
str(logist_binom_predict)
##  ppd [1:500, 1:19] 0.853 0.845 0.829 0.848 0.863 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:19] "1" "2" "3" "4" ...

\[ p(y^{rep}|y) = \int_{\Theta} p(y^{rep}|\theta) p(\theta|y)d\theta \]

A Geometric Model

\[ P(\text{success from dist } x) \\ = 2\Phi\Bigg(\frac{1}{\sigma} arcsin\Big(\frac{R-r}{x}\Big)\Bigg)-1 \]

P_success <- function(sigma, R, r, x){
  2 * pnorm(asin((R-r)/x) / sigma) - 1
  }
x <- seq(0,30,0.1)
qplot(x, P_success(1*pi/180, R, r, x),
      geom="line") + ylim(0,1)

Stan Model Part 1

data {
  int N;
  int<lower=0> tries[N];
  int<lower=0> successes[N];
  real<lower=0> dist[N];
  real R;
  real r;
}

parameters {
  real<lower = 0> sigma;
}

Stan Model Part 2

model {
  real p[N];
  
  for (n in 1:N) 
    p[n] = 2 * Phi(asin((R - r) / dist[n]) / sigma) - 1;
  
  successes ~ binomial(tries, p);
}

Vectorised:

model {
  real p[N];
  p = 2 * Phi(asin((R - r) / dist) / sigma) - 1;
  
  // A comment
  successes ~ binomial(tries, p);
}

Stan Model Part 3

generated quantities {
  vector[N] sucess_predictions;
  for (n in 1:N) {
    real p;
    p = 2 * Phi(asin((R - r) / dist[n]) / sigma) - 1;
    sucess_predictions[n] = binomial_rng(tries[n], p);
  }
}

Stan Data Types

int N;
real x;
vector[N] x;
simplex[K] x;
unit_vector[K] x;
ordered[K] x;
positive_ordered[K] x;

cov_matrix[K] x;
corr_matrix[K] x;

cholesky_factor_cov[K] x;
cholesky_factor_corr[K] x;

Prepare Data & Compile Model

str(data)
## List of 6
##  $ N        : num 19
##  $ tries    : num [1:19] 1443 694 455 353 272 ...
##  $ successes: num [1:19] 1346 577 337 208 149 ...
##  $ dist     : int [1:19] 2 3 4 5 6 7 8 9 10 11 ...
##  $ r        : num 0.07
##  $ R        : num 0.177
library(rstan)
#options(mc.cores = parallel::detectCores())
golf_stan <- stan_model(file = "golf.stan")
## In file included from filecd95fc1d82f.cpp:8:
## In file included from /Users/wells/R_libraries/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Users/wells/R_libraries/StanHeaders/include/stan/math.hpp:4:
## In file included from /Users/wells/R_libraries/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Users/wells/R_libraries/StanHeaders/include/stan/math/rev/core.hpp:12:
## In file included from /Users/wells/R_libraries/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5:
## In file included from /Users/wells/R_libraries/StanHeaders/include/stan/math/rev/core/var.hpp:7:
## In file included from /Users/wells/R_libraries/BH/include/boost/math/tools/config.hpp:13:
## In file included from /Users/wells/R_libraries/BH/include/boost/config.hpp:39:
## /Users/wells/R_libraries/BH/include/boost/config/compiler/clang.hpp:200:11: warning: 'BOOST_NO_CXX11_RVALUE_REFERENCES' macro redefined [-Wmacro-redefined]
## #  define BOOST_NO_CXX11_RVALUE_REFERENCES
##           ^
## <command line>:6:9: note: previous definition is here
## #define BOOST_NO_CXX11_RVALUE_REFERENCES 1
##         ^
## 1 warning generated.

Sample!

golf_fit <- sampling(golf_stan, data=data) 
## 
## SAMPLING FOR MODEL 'golf' NOW (CHAIN 1).
## 
## Gradient evaluation took 2e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.2 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.025384 seconds (Warm-up)
##                0.025904 seconds (Sampling)
##                0.051288 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'golf' NOW (CHAIN 2).
## 
## Gradient evaluation took 7e-06 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.025744 seconds (Warm-up)
##                0.026367 seconds (Sampling)
##                0.052111 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'golf' NOW (CHAIN 3).
## 
## Gradient evaluation took 8e-06 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.025561 seconds (Warm-up)
##                0.024421 seconds (Sampling)
##                0.049982 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'golf' NOW (CHAIN 4).
## 
## Gradient evaluation took 7e-06 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:    1 / 2000 [  0%]  (Warmup)
## Iteration:  200 / 2000 [ 10%]  (Warmup)
## Iteration:  400 / 2000 [ 20%]  (Warmup)
## Iteration:  600 / 2000 [ 30%]  (Warmup)
## Iteration:  800 / 2000 [ 40%]  (Warmup)
## Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Iteration: 2000 / 2000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.025476 seconds (Warm-up)
##                0.024774 seconds (Sampling)
##                0.05025 seconds (Total)

Sampling Options

golf_fit <- stan(file = "golf.stan",
                 model_name = "golf_putting",
                 data = data,
                 chains = 3,
                 iter = 1000,
                 warmup = 200,
                 thin = 2,
                 seed = 42,
                 cores = 2,
                 algorithm = "HMC")

golf_fit <- vb(golf_stan, data=data)
make golf
./golf sample data file=golf_data.R
./golf variational data file=golf_data.R

print(golf_fit)
## Inference for Stan model: golf.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##                            mean se_mean    sd     2.5%      25%      50%
## sigma                      0.03    0.00  0.00     0.03     0.03     0.03
## sucess_predictions[1]   1378.59    0.17  9.00  1360.00  1373.00  1379.00
## sucess_predictions[2]    568.41    0.19 11.08   547.00   561.00   568.00
## sucess_predictions[3]    311.62    0.18 10.37   291.00   305.00   312.00
## sucess_predictions[4]    203.93    0.16  9.76   185.00   197.00   204.00
## sucess_predictions[5]    135.01    0.14  8.34   119.00   129.00   135.00
## sucess_predictions[6]    111.08    0.13  8.16    95.00   106.00   111.00
## sucess_predictions[7]     92.38    0.13  7.63    77.00    87.00    92.00
## sucess_predictions[8]     74.92    0.11  6.98    62.00    70.00    75.00
## sucess_predictions[9]     62.37    0.11  6.70    50.00    58.00    62.00
## sucess_predictions[10]    67.47    0.12  6.98    54.00    63.00    67.00
## sucess_predictions[11]    52.98    0.10  6.31    41.00    49.00    53.00
## sucess_predictions[12]    46.54    0.10  5.99    35.00    43.00    47.00
## sucess_predictions[13]    39.26    0.09  5.50    29.00    36.00    39.00
## sucess_predictions[14]    35.16    0.08  5.28    25.00    32.00    35.00
## sucess_predictions[15]    39.88    0.09  5.61    29.00    36.00    40.00
## sucess_predictions[16]    36.41    0.09  5.40    26.00    33.00    36.00
## sucess_predictions[17]    33.74    0.08  5.28    24.00    30.00    34.00
## sucess_predictions[18]    24.74    0.07  4.53    16.00    22.00    25.00
## sucess_predictions[19]    24.14    0.07  4.50    16.00    21.00    24.00
## lp__                   -2926.77    0.02  0.72 -2928.80 -2926.92 -2926.48
##                             75%    97.5% n_eff Rhat
## sigma                      0.03     0.03  1719    1
## sucess_predictions[1]   1385.00  1396.00  2699    1
## sucess_predictions[2]    576.00   590.00  3270    1
## sucess_predictions[3]    319.00   331.03  3227    1
## sucess_predictions[4]    210.25   223.00  3691    1
## sucess_predictions[5]    141.00   151.00  3473    1
## sucess_predictions[6]    116.00   127.00  3790    1
## sucess_predictions[7]     97.00   107.00  3652    1
## sucess_predictions[8]     80.00    89.00  4000    1
## sucess_predictions[9]     67.00    76.00  4000    1
## sucess_predictions[10]    72.00    81.00  3678    1
## sucess_predictions[11]    57.00    65.00  3957    1
## sucess_predictions[12]    50.00    59.00  3953    1
## sucess_predictions[13]    43.00    50.00  4000    1
## sucess_predictions[14]    39.00    46.00  4000    1
## sucess_predictions[15]    44.00    51.00  4000    1
## sucess_predictions[16]    40.00    47.00  3681    1
## sucess_predictions[17]    37.00    44.00  4000    1
## sucess_predictions[18]    28.00    34.00  4000    1
## sucess_predictions[19]    27.00    33.00  3944    1
## lp__                   -2926.31 -2926.27  1671    1
## 
## Samples were drawn using NUTS(diag_e) at Sun Feb  4 00:56:23 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

Diagnostics

stan_trace(golf_fit, pars = "sigma")

stan_ac(golf_fit, pars = 'sigma')

Extract Samples

str(as.matrix(golf_fit))
##  num [1:4000, 1:21] 0.0267 0.0266 0.0265 0.0266 0.0267 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ iterations: NULL
##   ..$ parameters: chr [1:21] "sigma" "sucess_predictions[1]" "sucess_predictions[2]" "sucess_predictions[3]" ...
str(as.array(golf_fit))
##  num [1:1000, 1:4, 1:21] 0.0267 0.0266 0.0265 0.0266 0.0267 ...
##  - attr(*, "dimnames")=List of 3
##   ..$ iterations: NULL
##   ..$ chains    : chr [1:4] "chain:1" "chain:2" "chain:3" "chain:4"
##   ..$ parameters: chr [1:21] "sigma" "sucess_predictions[1]" "sucess_predictions[2]" "sucess_predictions[3]" ...

str(rstan::extract(golf_fit))
## List of 3
##  $ sigma             : num [1:4000(1d)] 0.027 0.027 0.0261 0.0272 0.0273 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL
##  $ sucess_predictions: num [1:4000, 1:19] 1364 1382 1384 1372 1363 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ iterations: NULL
##   .. ..$           : NULL
##  $ lp__              : num [1:4000(1d)] -2927 -2927 -2927 -2927 -2927 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL

Check Fit

Posterior Predictive

One liner PPCs with library(bayesplot)

ppc_violin_grouped(y=golf$p,
  yrep=logist_binom_predict[1:500,],
  group = golf$dist, y_size = 2, y_draw = "both") +
  theme(legend.position = "bottom")

ppc_violin_grouped(y=golf$p,
  yrep=golf_predictive[1:500,],
  group = golf$dist, y_size = 2, y_draw = "both") +
  theme(legend.position = "bottom")

ppc_stat_grouped(golf$p,
  logist_binom_predict[1:500,],
  group = golf$dist, stat = "median") +
  theme(legend.position = "bottom")

ppc_stat_grouped(golf$p,
  golf_predictive[1:500,],
  group = golf$dist, stat = "median") +
  theme(legend.position = "bottom")

mcmc_areas(as.matrix(golf_fit)[,c(11,12,13), drop=FALSE])

ppc_dens_overlay(golf$p,
                 golf_predictive[1:500,])

ppc_dens_overlay(golf$p,
                 logist_binom_predict[1:500,])

Shinystan

library(shinystan)
launch_shinystan(golf_fit)

Alternative 'Sampling' Statements

beta ~ normal(0, 1);
target += normal_lpdf(beta | 0, 1);
target += -0.5 * beta^2;
target += -0.5 * beta’ * beta;
target += -0.5 * dot_self(beta);

Stan Code Linear Regression

data {
  int           N ; # integer, number of observations
  int           K ; # integer, number of columns in model matrix
  matrix[N,K]   X ; # N by K model matrix
  vector[N]     y ; # vector of N observations
}

parameters {
  real<lower=0> sigma ; # real number > 0, standard deviation
  vector[K]     beta ;  # K-vector of regression coefficients
}

model {
  beta ~ normal(0, 5) ;       # prior for betas
  sigma ~ cauchy(0, 2.5) ;    # prior for sigma
  y ~ normal(X*beta, sigma) ; # vectorized likelihood
}

generated quantities {
  vector[N] y_rep ; # vector of same length as the data y
  for (n in 1:N) 
    y_rep[n] <- normal_rng(X[n]*beta, sigma) ;
}

Refrences